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Outline 


Q  Plasma 

Q  Multi-Fluid  Plasma  Model 
o  Advantages 
o  Limits  of  the  model 

Q  Numerics:  Blended  Finite  Element  Method 
o  Discontinuous  Galerkin  for  ion/neutrals 
o  Continuous  Galerkin  for  electrons/fields 
o  Initial  tests 
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What  are  plasmas? 


Plasma  is  a  quasineutral  gas  of  charged  and  neutral  particles  which  exhibits 

collective  behavior. 


Credit:  particleincell.com 


o  “99%  of  matter  in  the  universe  is  in  the  state 
of  plasma" 
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Time/Spatial  scales. 


□ 


Sean  Miller,  PhD  dissertation,  University  of  Washington  (2016) 
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Velocity 


There  are  multiple  plasma  models. 


o  3-Dimensions  +  3-Velocities 

o  Evolve  the  particles  position 
and  velocity 

o  e.g.  Particle-In-Cell  models 


a  Ensemble  average  of  particles 
distribution,  fs(x,x,t) 

o  Evolve  the  distribution 
function 

o  e.g.  Vlasov-Maxwell  models 
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Boltzmann  equation  e voles  fs 


SP 


o  The  Boltzmann  eqn: 


dfs  dfs  qs  .  <9/v  dfs 

dt  dx  msK  J  d\  dt 


o  Take  the  0th ,  kst,  2/W  moments  of  the  Boltzmann  Eqn. 


v”f  dy+m>\ 


n+1 


ax 


d\  +  qs  J  \n  (E  +  v  X  B)  •  =  ms  J  \n 


% 

dt 


d\ 


o  Each  moment  of  the  Boltzmann  eqn  gives  an  equation  for  the  moment 
variable,  and  introduces  the  next  higher  moment  variable 

o  This  process  can  go  on  indefinitely 
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Boltzmann  equation  e voles  fs 


dps 

dt 


+  V  •  (ps us) 


dps 

dt 


dps  u* 
dt 


+  V  •  (ftUslli  +  Ps  I  +  IIS) 


^(E  +  UjxB)-VRJiW  + 
ms 


dps  us 
dt 


r 


des 

dt 


+  V  •  (((es  +  Ps)  I  +  n s)  •  u,  +  hs) 


— u,  •  E  +  Qs  s* 

m.  z ' 


des 

dt 


r 


o  System  is  truncated  by  relating  higher  moment  variables  to  the  lower 
ones 

o  The  fluids  are  coupled  to  each  other  and  to  the  electromagnetic  fields 
through  Maxwell’s  equations  and  interaction  source  terms. 


Sousa  (ERC/AFRL) 


Distribution  A:  Approved  for  public  release;  unlimited  distribution  cieara 


No.  17211 


7/31 


Plasma  models  range  of  applicability. 


KJ 


I"gi0  (n^) 


Q 


Sean  Miller,  PhD  dissertation,  University  of  Washington  (2016) 
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Advantages  of  the  Model 


LTE,  velocity  moments  ,  ,  ea— M3,  me—t 0  m  ,TTTX 

Kinetic  - - - >  MFPM  — — - - — >  MHD 


Ideal  MHD  model  is  valid  when: 


o  High  collisionality,  tu/t  <C  1 
o  Small  Larmor  radius,  rLi/L  <C  1 

/  \l/2  9 

*  Low  Resistivity,  (^)  (?)  %«1 


Multi-Fluid  Plasma  Model 


9  Less  computationally  expensive  than  kinetic  models 
o  Multi-fluid  effects  become  relevant  at  small  spacial  and  temporal  scales 
o  Finite  electron  mass  and  speed-of-light  effects  are  included 
o  There  is  charge  separation  is  modeled 
o  Displacement  current  effects  are  resolved  in  the  MFPM 
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The  MFPM  has  dispersive  sources. 


dQ  dV  _ 

dt  dx 


The  source  Jacobian  has  imaginary  eigenvalues 

o  Q 

The  equation  system  has  dispersive  sources 

The  dispersion  is  physical  (may  be  difficult  to  distinguish  from 
numerical  dispersion) 


9  This  dispersion  is  due  to  plasma  waves  that  result  from  ion  and  electron 
plasma  interactions  with  electromagnetic  fields 

9  An  ideal  numerical  method  for  the  MFPM  should: 

9  be  high-order  accurate  9  couple  the  flux  and  the  sources 

9  capture  shocks  o  not  impose  strict  time-step 

j|  Srinivasan  et  al,  CCP  10  (2011) 
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BFEM  SIMULTANEOUSLY  USES  CG  AND  DG. 


9  Solution  to  the  electron  and  EM  fields  is  smooth  and  does  not  shock 


9  Continuous  Galerkin 


9  Discontinuous  Galerkin 


9  Electron  fluid  and  EM  fields 


9  Multiple  ion  and  neutral  fluids 


Q  =  ^q/v« 

i 


q  =  J2cm 

i 
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Implicit  Continuous  Galerkin 


o  For  this  implementation  the  balance  law  form  is  cast  as 

dQ  d ^  dQ  9 

^  +  W^  =  s  +  KVQrf 

o  Lagrange  polynomials  are  used  for  basis  functions,  vr 

dr-  ~r(Q)  =  f  vrSdV  -  f  9F 
Jn  ai  Jn  Jn 

o  # -method  time  integration 

=  M^"+1  “  Q" 


/ 

Jn 


vr—^—dV  =  Cr(Q)  =  [  vrS dV  —  [  dV  +  k  f  vrV2Q ddV 

Jn 


dQ  dx 


nor) 


dt 


-  0Cr(Qn+l)  -  (1  -  0)Cr(Qn)  =  0 


o  6  =  0.5  is  used  for  2nd  order  accuracy 

?(Qn)  =  ^^,  Y  (Q")AQ  =  -K(Qn),  Q”+1  =  Q”  +  AQ 

jj  Reddy,  An  Introduction  to  the  Finite  Element  Method  (2006) 
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Runge-Kutta  Discontinuous  Galerkin 


9Q  ,  dV  _ 

dt  dx 

o  Legendre  polynomials  are  used  for  basis  functions,  vp 
o  The  hyperbolic  equation  is  multiplied  by  the  basis  function, 

[  vp^dV  =  Cp( Q)  =  [  vpSdV  -  (f  Vp^F  •  dA  +  [  V  •  X7vpdV 

Jn  dt  Jq  Jqq  Jq 

©  Explicit  Runge-Kutta  time  integration 
©  CFL  =  cAt/Ax  <1/ (2p  —  1),  /7  is  the  polynomial  order 

Q*  =  Q"  +  AtCp(Qn), 

Qn+i  =  Iq*  +  Iq«  +  Ia^(q*). 

j|  Loverich  et  al,  CCP  9  (2006) 
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Convergence  of  the  BFEM. 


9  Simulations  at  fixed  time- step _ 9  Simulations  at  fixed  CFL=1 
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ID  Soliton  Problem 


a  ID  soliton  is  a  two-fluid  plasma 
problem 

o  The  solution  is  smooth,  therefore 
artificial  dissipation  can  be  small 

o  The  simulation  uses  512 
second-order  elements 
o  Bz  =  1.0,  Te  =  Ti  =  0.01, 

U/  =  u,  =  0 

9  ne  =  m  —  1.0  +  e_10(x_6)2 

j]  Baboolal,  Math,  and  Comp.  Sim.  55  (2001) 
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DG  AND  BFEM  COMPARISON  WITH  A  SOLUTION 


—  =  1836,  —  =  1000V2,  FV  5000  cells 

csi 


9  DG  Solution  is  very  dispersive 

Hi  Hakim  et  al,  JCP  219  (2006) 


®  BFEM  is  less  dissipative  than  the 
converged  solution 


Sousa  (ERC/AFRL) 


Distribution  A:  Approved  for  public  release;  unlimited  distribution  clearance  no.  17211 


16/31 


BFEM  Computational  Cost  Savings 


case 

mj/me 

C/Csi 

DG  time(s) 

BFEM  time  (s) 

BFEM  cost  over  DG 

1 

25 

10/72 

0.32 

37.7 

+11681% 

2 

100 

10/72 

1.28 

37.7 

+2845% 

3 

500 

10/72 

6.82 

37.7 

+452.8% 

4 

1000 

10/72 

12.4 

38.2 

+208.1% 

5 

1836 

10/72 

23.5 

40.4 

+71.91% 

6 

3672 

10/72 

47.2 

39.2 

-16.95% 

7 

3672 

100/72 

520 

265 

-49.04% 

8 

3672 

1000/72 

5274 

2735 

-48.14% 

o  Per  time-step  explicit  DG  is  faster  than  BFEM,  but  it  requires  many  more 
time-steps 

o  BFEM  is  more  efficient  only  when  time-step  are  considerably  larger  than 
explicit  DG 
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L2-norm 


BFEM  Accuracy 


o  The  BFEM  seems  to  be  less  accurate  than  the  DG  implementation 
(~  50%) 

o  When  the  mass  ratio  is  one,  the  two  methods  have  the  same  level  of 
accuracy 

o  The  discrepancy  is  due  to  the  fact  that  the  semi-implicit  BFEM  does  not 
resolve  the  plasma  frequency  in  this  problem 
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Electromagnetic  Plasma  Shock  Problem 


a  Fast  rarefaction  wave  (FR),  a 
slow  compound  wave  (SC),  a 
contact  discontinuity  (CD),  a 
slow  shock  (SS),  and  another 
fast  rarefaction  wave  (FR) 

a  The  problem  exhibits  limits 
of  MHD  and  multi-fluid 
behavior  by  changing  the 
Larmor  radius,  r l 
o  MHD:  rL  — 0 
o  Multi-fluid:  vL 

j|  Brio  and  Wu,  JCP  75  (1988) 
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Shock  in  Density  but  Smooth  Fields. 


o  t=0.05/cca-,  c/c^-110,  mz7nig=1836 


mass  density  Magnetic  field  (y-comp.) 

o  The  main  features  of  the  problem  are  captured  by  all  three  methods 

o  BFEM  does  not  properly  resolve  the  fast  electromagnetic  waves  which 
require  accurately  resolving  the  electron  dynamics 
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Maximum  BFEM  time-step 


A tmax  —  min 


Ax  Ax  Ax  0.1  0.1  0.1  0.1 


Ur, 


Uci  Upe  UJpi 


9  A tmax  corresponds  to  the 
maximum  value  allowed  for 
explicit  methods  based  on  the 
CFL  condition 

9  At  =  42.9A tmax  is  the 
maximum  time  step  allowed 
by  the  BFEM  due  to  ion 
dynamics 
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Effects  of  Artificial  Dissipation 


9  Varying  the  artificial 
dissipation  on  the  electron 
fluid,  Ke 

9  Wave-like  behavior  of  the 
problem  is  better  resolved 

9  Amplitude  of  the  compound 
wave  increases 

9  Right  fast  rarefaction  wave  is 
not  visible 
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Effects  of  Artificial  Dissipation 


®  Varying  the  artificial 

dissipation  on  the  EM-field, 

^EM 

9  There  is  better  agreement 
with  the  DG  solution 

o  This  reinforces  the  point  that 
the  wave-like  behavior  arises 
from  the  interaction  of  the 
electron  fluid  with  the 
electromagnetic  fields 
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Summary 


o  The  blended  finite  element  method  (BFEM)  is  presented 
o  DG  spatial  discretization  with  explicit  Runge-Kutta  (i+,  n) 
o  CG  spatial  discretization  with  implicit  Crank-Nicolson  (e-,  fileds) 
o  DG  captures  shocks  and  discontinuities 
o  CG  is  efficient  and  robust  for  smooth  solutions 

o  Physics-based  decomposition  of  the  algorithm  yields  numerical  solutions 
that  resolve  the  desired  timescales 

e  DG  method  takes  less  computational  time  to  advance  the  solution  by  one 
time-step,  however  At  is  much  smaller  than  that  of  the  BFEM 

o  Computational  cost  savings  using  the  BFEM  will  only  occur  for 
relatively  large  implicit  time-steps  compared  to  explicit  time-steps 

%  Sousa  and  Shumlak,  JCP  326  (2016)  56-75 

Thank  you. 
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EXTRAS 
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Physical  description  of  a  fluid 


o  Modeling  each  particle  velocity  and  position  is  not  practical. 

©  Instead  an  average  is  performed  to  give  a  statistical  description. 

o  Calculate  the  number  of  particles  per  unit  volume  having  approximately 
the  velocity  v  near  the  position  x  and  at  time  t ,  distribution  function 

/(v,x,  t) 


Ps  =  Ms 
psus  =  ms 


J  fs(v)d\ 
j  \fs(\)d\ 

=  P.v  =  ms  J  ww fs(v)dv,  ps  =  J  w2fs(\)d\ 
J  WWW fs(v)d\,  hv  =  ]-ms  J  w2wfs(\)d\ 


Ho  =  m 


w  =  v  -  Uo 
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Field  Reversed  Configuration  Thruster  o 


1  -  Formation  of  high  density  FRC 

2  -  Acceleration  of  the  FRC  by 
Lorentz  force 

3  -  FRC  expands,  converting 
thermal  energy  to  directed  energy 


Advantages: 

•  Electrodeless  and  the  plasmoid  propellant  is  magnetically  isolated  from  the  walls 

•  Propellant  is  completely  uncoupled  from  the  driving  and  confining  fields 

•  High  plasma  temperatures  and  densities  significantly  reduce  ionization  losses 
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FRC  FORMATION  STUDIES. 


9  The  formation  of  FRC  using  a  rotating  magnetic  field  (RMF) 

O  LOci  <U  <iL)ce 
9  Vei  <  CCce 


9  Flux  Conserver  BC 

b2Ba{  0)  1 


Bz(a,t)  = 


b2  —  a2  7 r(b2  —  a2) 


fa  rdr  f2 

Jo  Jo 


d6Bz(r,0,t) 


o  RMF  BC 

Irmf  =  It  cos(a Jt  +  (/>) 
I,  =  I0(l-e~,/T) 

®  4>A,B  =  0,  4>e,F  =  7T 
®  4>c,d  =  tt/2, 

=  -'ll  1 1 

9  Bias  Field  is  constant  at 
t=0s 
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RMF  Antenna  Modeling 


©  Initialization:  B^WG,  BZ=50G,  cc=5MHz 
o  Solve  Maxwell’s  eqns.  with  divergence  constraints 
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FRC  Formation  Simulations 


o  Initialization:  T;=Tg=30eV,  n/=ng=1019m-3,  m//me=1836 
o  Solve  Multi-Fluid  eqns. 
o  Bz  is  plotted 
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(Near)  Future  Work 


Neutral-Plasma  Model 

o  Interaction  of  multi-fluid  plasma 
with  gas  dynamics  neutral  fluid 
9  electron-impact  ionization 

e  +  n  — >  +  7,e  —  (pion 

9  radiative  recombination 
e~  +  i+  — >•  n  +  hv 
9  resonant  charge  exchange 
z+  +  n  — >>  n  +  /+ 

j)  Meier  and  Shumlak,  PoP  19  072508  (2012) 
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Collisional-Radiative  Model 


o  Excitation/De-excitation  rates 


Q 


Le  and  Cambier,  PoP  22  093512  (2015) 


